gc()

setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
setwd("../../")

## Load/install packages
if (!require("pacman")) install.packages("pacman")
pacman::p_load(tidyverse,
               parallel,
               ggplot2,
               data.table,
               viridis,
               miceadds,
               haven,
               cdlTools,
               rgdal,
               sf,
               tidycensus)

savefig <- function(path,width,height) {
  
  lapply(c("png","pdf"), function(ftype) 
    ggplot2::ggsave(paste0(path,".",ftype),
                    device=ftype,width=width,height=height))
  
}

western_states <- c("AZ","CA","CO","ID","MT","NM","NV","OR","WA","WY","UT")
states <- read_sf("data/geo/western_states/western_states.shp")

coords <- states %>% filter(STUSPS == "CO") %>% st_transform(st_crs(4269)) %>% st_union() %>% st_centroid() %>% st_coordinates()
crs <- sprintf("+proj=aea +lat_1=29.5 +lat_2=45.5 +lat_0=%s +lon_0=%s +x_0=0 +y_0=0 +datum=NAD83 +units=m +no_defs", coords[[2]], coords[[1]])
  
states <- st_transform(states, st_crs(crs))

get_fires <- function(sample) {
  
  ignition_ids <- duration_data %>% filter(firenum %in% sample)
  ignition_ids <- ignition_ids$ignitionid
  
  df <- read_dta("data/IgnTableWest.dta") %>%
    dplyr::select(c("IgnitionID","STATE","OWNER_DESCR","FIRE_YEAR","FIRE_SIZE","DISCOVERY_DATE","STAT_CAUSE_DESCR",
                    "LATITUDE","LONGITUDE")) %>%
    filter(IgnitionID %in% ignition_ids)
  
  df <- st_as_sf(df, coords = c("LONGITUDE","LATITUDE"), remove = FALSE, crs = 4269) %>%
          st_transform(st_crs(crs))
  
  return(df)
  
}

duration_data <- read_dta("data/20-by-1km_15deg/data_centroid.dta") %>% 
  group_by(firenum,ignitionid) %>% 
  summarize(num = sum(1))

sample.df <- read_dta("data/20-by-1km_15deg/sample_census.dta") 
sample_census <- sample.df %>% filter(sample == 1)
sample_census <- unique(sample_census$firenum)
df_census <- get_fires(sample_census)

sample.df <- read_dta("data/20-by-1km_15deg/sample_assessors.dta") 
sample_assessors <- sample.df %>% filter(sample == 1)
sample_assessors <- unique(sample_assessors$firenum)
df_assessors <- get_fires(sample_assessors)

fires <- rbind(df_census %>% mutate(sample = "Census sample"),
               df_assessors %>% mutate(sample = "2012-2015 assessors' data sample")) %>%
  mutate(sample = fct_relevel(sample, "Census sample", "2012-2015 assessors' data sample"))


ggplot(data = states) +
  geom_sf() + 
  geom_sf(data = fires, aes(color = sample, shape = sample)) + 
  labs(color = "Ignition locations") +
  scale_color_manual("Ignition locations", labels = c("Census sample","2012-2015 assessors' data sample"), values = c("black","red")) +
  scale_shape_manual("Ignition locations", labels = c("Census sample","2012-2015 assessors' data sample"), values = c(16,17)) +
  theme(plot.background = element_rect(fill = "transparent",color = NA),
        panel.border = element_blank(),
        panel.background =element_rect(fill = "transparent",color = NA),
        panel.grid = element_blank(),
        plot.margin = unit(c(0,0.25,0,0), "cm"),
        plot.subtitle = element_text(hjust = 0.5),
        axis.text = element_blank(),
        axis.ticks = element_blank(),
        axis.title = element_blank(),
        legend.direction = "vertical")

dst = "Results/western_fire_sample_map_vR"
savefig(dst, width = 8*1.25, height = 5*1.25)
